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In general neutron stars in binaries are spinning. Due to the existence of millisecond pulsars we 
know that these spins can be substantial. We argue that spins with periods on the order a few dozen 
milliseconds could influence the late inspiral and merger dynamics. Thus numerical simulations of 
the last few orbits and the merger should start from initial conditions that allow for arbitrary spins. 
We discuss quasi-equilibrium approximations one can make in the construction of binary neutron 
star initial data with spins. Using these approximations we are able to derive two new matter 
equations. As in the case of irrotational neutron star binaries one of these equations is algebraic 
and the other elliptic. If these new matter equations are solved together with the equations for 
the metric variables following the Wilson-Mathews or conformal thin sandwich approach one can 
construct neutron star initial data. The spin of each star is described by a rotational velocity that 
can be chosen freely so that one can create stars in arbitrary rotation states. Our new matter 
equations reduce to the well known limits of both corotating and irrotational neutron star binaries. 



PACS numbers: 04.20.Ex, 04.30.Db, 97.60.Jd, 97.80.Fk 



I. INTRODUCTION 

Several gravitational wave detectors such as LIGO [l|, 
0] , Virgo ^3, 4] or GEO @ have been operating over the 
last few years, while several others are in the planning 
or construction phase [y]. One of the most promising 
sources for these detectors are the inspirals and mergers 
of binary neutron stars. In order to make predictions 
about the last few orbits and the merger of such sys- 
tems, fully non-linear numerical simulations of the Ein- 
stein Equations are required. To start such simulations 
we need initial data that describe the binary a few orbits 
before merger. The emission of gravitational waves tends 
to circularize the orbits |7|, |8|. Thus, during the inspiral, 
we expect the two neutron stars to be in quasi-circular 
orbits around each other with a radius that shrinks on 
a timescale much larger than the orbital timescale. This 
means that the initial data should have an approximate 
helical Killing vector f^. To incorporate these ideas we 
will use the Wilson-Mathews approach [3, Il3] ) which is 
also known as conformal thin sandwich formalism [llj, 
for the metric variables. The Wilson-Mathews approach 
has already been successfully used by several groups to- 
gether with matter equations describing the neutron stars 
in either corotating |13 - [l6| or irrotational |17H25| states. 
There have also been attempts to include intermediate 
rotation states [2l], [2g]. However, as we will discuss in 
more detail later, these approaches have certain draw- 
backs, because they do not correctly solve the Euler equa- 
tion for the fluid. Thus, so far there is no canonical for- 
malism to describe neutron star binaries with arbitrary 
spins. As pointed out by Bildsten and Cutler [23|, the two 
neutron stars cannot be tidal locked, because the viscos- 
ity of neutron star matter is too low. Hence barring other 
effects like magnetic dipole radiation the spin of each star 
remains approximately constant. This means that initial 
data sequences of corotating configurations for different 
separations cannot be used to approximate the inspiral 
of two neutron stars. On the other hand, sequences of 



irrotational configurations can be used to approximate 
the inspiral of two neutron stars without spin. This fact 
explains why irrotational initial data are far more popu- 
lar today. Nevertheless, astrophysical neutron stars will 
have a non-zero spin. Therefore a corotating configura- 
tion at some particular separation does have its place as 
a possible initial configuration with spin. It will just not 
remain corotating during the subsequent time evolution. 
Of course real neutron stars will likely have spins that 
have periods different from the orbital period, and the 
spin direction may not be aligned with the orbital angu- 
lar momentum. Thus it would be highly desirable to have 
a formalism that can be used to generate initial data for 
arbitrary initial spins. 

In order to judge how important spins might be let us 
discuss a few order of magnitude estimates. A typical 
neutron star has a mass of about 1.4 solar masses (Mq) 
and a radius on the order of I5km. From Kepler's law 
the orbital period 
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is on the order of a few milliseconds during the last or- 
bit before merger where the separation d ~ 50km. Thus 
systems with spin periods that are much larger than Pq 
should be treatable as approximately irrotational, while 
systems with spin periods of a few milliseconds (such as 
millisecond pulsars) cannot be regarded as irrotational. 
Another way of judging how important spins could be 
during the evolution is to look at the dimensionless spin 
magnitude. If we assume that the spin 5* of a neutron 
star with mass m and radius R is related to its spin pe- 
riod P by S" = /(27r/P) with / ~ mR"^ we find that the 
dimensionless spin has a magnitude of 
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Thus millisecond pulsars have a dimensionless spin of or- 
der one. As in the case of binary black holes |28l - [33| . 



spins of this magnitudes could have a significant influ- 
ence on the merger dynamics. This means that neutron 
stars with spin periods of a few dozen milhseconds or 
less should not be considered irrotational. One could of 
course imagine that the neutron stars spin down before 
they enter the strongly relativistic regime of the last few 
orbits before merger that is usually considered in numer- 
ical relativity simulations. In order to address this ques- 
tion let us look at the famous double pulsar PSR J0737- 
3039 which is the only neutron star binary where both 
spin periods and spin down rates are known [33] • Star 
A has mass ttia = 1.34Mq and spin period Pa = 23ms, 
while star B has tob = 1.25M0 and Pb = 2.8s. The 
orbital period is Po — 2.4h. From these numbers one 
derives that the system will merge in about 85My due 
to the emission of gravitational waves. Both stars are 
currently spinning down at a rate of Pa = 1-7 x 10^^® 
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34| . If one assumes that this spin 



down is due to magnetic dipole radiation and defines the 
characteristic ages given by ta = Pa/{2Pa) = 210My 
and tb — PbU'^Pb) = 50My, one finds that the spin 
period of each star obeys [35| 



PA/B{t)^PA/B{0)Jl 
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where the time i = is the time today. From this it is 
clear that the periods at merger (at t — 85My) will be 
PAit) ^ 27ms and PBit) = 4.6s. Thus star A will not 
spin down enough to be well approximated by an irro- 
tational configuration by this time. This example shows 
that neutron stars in binary systems can have apprecia- 
ble spins a few orbits before merger. Of course since 
only about ten binary neutron stars have been observed 
so far [36[ it is not clear yet how common binary neutron 
stars with high spins are. However, since there are nu- 
merous millisecond pulsars it seems reasonable to expect 
that neutron stars in binaries can also have millisecond 
spin periods. Hence the widely held belief that only ir- 
rotational configurations are realistic is not necessarily 
correct. It is thus necessary to develop initial data for 
binary neutron stars with arbitrary spins. In the next 
sections we will describe what approximation one can 
make to derive a formalism that allows for this possibil- 
ity. We will see that our new equations reduce to well 
known accepted results in both the corotating and irro- 
tational cases. 

Throughout we will use units where G = c = 1 . Latin 
indices such as i run from 1 to 3 and denote spatial in- 
dices, while Greek indices such as /x run from to 3 and 
denote spacetime indices. The paper is organized as fol- 
lows. Sec. [H] lists the General Relativistic equations that 
govern binary neutron stars described by perfect fluids. 
We use three approximate quasi-equilibrium conditions 
to simplify these equations. We find two new matter 
equations that allow as to set up binary neutron stars 
with arbitrary spins. In Sec. IIIII we consider the New- 
tonian limit of our new equations. We conclude with a 



discussion of our method in Sec. IIVI In the appendix 
we discuss our quasi-equilibrium conditions for a simple 
case. 



II. BINARY NEUTRON STARS WITH 
ARBITRARY ROTATION STATES 

In this section we describe the equations governing bi- 
nary neutron stars in arbitrary rotation states in Gen- 
eral Relativity. The equations for the metric and mat- 
ter variables discussed in subsections III Al IIIBl III CI and 
IIIDI are well known. Our new results concerning quasi- 
equilibrium conditions for neutron stars with arbitrary 
rotation states are presented in subsections III El and III Fl 



A. ADM decomposition of Einstein's equations 

We use the Arnowitt-Deser-Misner (ADM) decompo- 
sition of Einstein's equations (see e.g. (37|) and introduce 
the 3-metric 



7mi> = 5^"^ + '^m""- 



(4) 



Here g^jy is the spacetime metric and n^ is the unit nor- 
mal to the t ~ const hypersurface. The line element is 
then 

ds^ = -a^dt^ + 7y {dx' + I3'dt){dx' + pHt), (5) 

where the lapse a and shift /3* are related to n^ via 



n^ = (1/a, -/?Va) 



(-a, 0,0,0) (6) 



The extrinsic curvature is defined by 



(7) 



With these definitions Einstein's equations split into the 
evolution equations 






iirSij + AttJij {S - p) 



(8) 



and the Hamiltonian and momentum constraint equa- 
tions 






(9) 



Here Rij and R are the Ricci tensor and scalar computed 
from 7ij, Di is the derivative operator compatible with 
jij and all indices here are raised and lowered with the 
3-metric -fij. The source terms p, f, Sij and S = "f^-' Sij 
are projections of the stress-energy tensor T^^ given by 



p = T^^n'^n" 



(10) 



and correspond to the energy density, flux and stress- 
tensor. 



B. Matter equations 

We assume that the matter in both stars is a perfect 
fluid with a stress-energy tensor 

T^'' ^[p^{l + e) + P]u^'u" + Pg^'r (11) 

Here po is the mass density (which is proportional the 
number density of baryons), P is the pressure, e is the 
internal energy density divided by po a-nd u^ is the 4- 
velocity of the fluid. The matter variables in Eg. lfTO)) are 
then 

p = a^[po{l + () + P]u^vP-P 

f = a[po{l + e) + P]u\°{uyu° + l3') 

S'^ = [po(l + e) + P]w%0(wVuO + /3*)(uV"° + /3'') 

+PY' (12) 

From V^T^" = we obtain the relativistic Euler equa- 
tion 

[po(l + e) + P]u"V,u^ = ^{g^" + uf'u'')V,P, (13) 

which together with the continuity equation 

V.ipou") = (14) 

governs the fluid. 

In order simplify the problem we assume that inter- 
nal energy e is a function of po alone (which implies a 
temperature of zero), and use a polytropic equation of 
state 

P = .pJ+^/". (15) 

We also introduce the specific enthalpy 

h = l + e + P/po. (16) 

Changes in h at zero temperature obey 

dh = dP/po. (17) 



Using Eqs. p6)) and (|T7|) we can rewrite the Eulcr 
Eq. nil) as 



u'^V^u^ -H V^/i = 0, 



where 



(18) 



(19) 



It is often convenient to introduce the dimensionless 
ratio 



Q = P/Po, 
which we can use to write 

h = {n + l)q + 1 

Po - ^-"g" 

P = K-"g"+i 
e = nq. 



(20) 



(21) 



C. Decomposition of 3-nietric and extrinsic 
curvature 

As in [3, [iO| the 3-metric jij is decomposed into a 
conformal factor ^ and a conformal metric 7^ such that 



Itj 



Vlij- 



(22) 



The extrinsic curvature is split into its trace K and its 
tracefree part Aij by writing it as 



i^ij — Aij + —'jijK 



(23) 



D. 



Quasi-equilibrium assumptions for the metric 
variables 



We now make some additional simplifying assump- 
tions. First we assume that the binary is in an approxi- 
mately circular orbit and that the spins of each star re- 
main approximately constant. As in the case of binary 
black holes (see e.g. [3l,|33|) this implies the existence of 
an approximate helical Killing vector ^'' with =f ^5^1/ ~ 0. 
In order to clarify the meaning of the approximate sign 
we now briefly discuss two cases. 

If both spins are parallel to the orbital angular mo- 
mentum we have ££^g^v = 0{Po/Tins), where we assume 
the inspiral timescale Tins to be much longer than the 
orbital timescale Po- I.e. in a corotating coordinate sys- 
tem all metric time derivatives are of order 0{Po/Tins) 
and thus small. For arbitrary spins the situation becomes 
more complicated. We can again use corotating coordi- 
nates, but in this coordinate system the spin vectors will 
be precessing on an orbital timescale Pg- This means 
there are matter currents that change on a timescale Po, 
while the matter distribution itself only changes on the 
inspiral timescale Tins- In this case it is useful to con- 
sider gravity to be made up of gravitoelectric and grav- 
itomagnetic fields [33, |40[- The gravitoelectric parts of 
the metric are sourced by the matter distribution and 
thus change only on the timescale Tins, while the grav- 
itomagnetic parts of the metric are sourced by matter 
currents and thus change on the shorter timescale Po- 
However, the gravitomagnetic parts are smaller than the 
gravitoelectric parts by 0{v/c) [33, |40|- Thus we now 
have £^g^i, = 0{v/c) « 0, where we assume that the 
orbital velocity v is smaller than the speed of light. 

An approximate helical Killing vector with £^g^i, ~ 
implies that 



£^11 



£fK kO, 



(24) 



which is what we will need to assume here for the metric 
variables. In a corotating coordinate system where the 
time evolution vector lies along £^^, the time derivatives 
of these metric variables are then equal to zero. From 
dtjij = it follows that 



A'^ 



1 



2iP^a 



[Lpy- 



(25) 



where 

{Lpy^ = &p^ + &p' - lDk/3\ (26) 

and D^. is the derivative operator compatible with 7,^. 
The assumption dtK = together with the evolution 
equation of K (derived from Eq. ([5])) implies 

+47ra(S'-3p). (27) 



E. Quasi-equilibrium assumptions for the matter 
variables 



In an inertial frame (i.e. a frame with limr^.oo P^ = 0) 
the approximate helical Killing vector has the compo- 
nents 



e 



[\,-n\x^ 



^Cm\ ; 



VL\x^ 



''CMJ; 



0) 



(28) 



Here 



'CM 



denotes the center of mass position of the sys- 



tem (which can be obtained from surface integrals at 
infinity e.g. Eq. (20.11) in (s^l), and D, is the orbital 
angular velocity, which we have chosen to lie along the 
a;"^-direction. Following Shibata [41| we decompose the 
fluid velocity u'^ into a piece along ^^ and a spatial vec- 
tor V^ and write 



u^' = u° [e + V^") 



(29) 



where vP — —u^n^/ a. 

In terms of ^'^ and V^^ the fluid equations p^ and ([1 
can be recast as 

A (poauV*) -f a [£^{pou'') + p^u"^ g^''' £ ^g ^,] = (30) 

and 

A (^ +^''^UkV'^^+V'' (ofu, - Dfkk)W^£^u, = 0, 

(31) 
where 



(3)m» = f^u\ 



(32) 



When one constructs neutron star initial data for coro- 
tating or irrotational configurations one usually assumes 
that the Lie derivatives of all matter variables with re- 
spect to ^'^ vanish [ij, |4l|, [IJ. However, for arbitrary 
spins this may not be the best approximation, since the 
portion of the fluid velocity responsible for the star's spin 
is not constant along ^^ if the spin remains constant while 
the stars orbit around each other. So we should not as- 
sume that £i^u^ vanishes. Rather we will split u^ into an 
irrotational and a rotational part and assume that only 
the Lie derivative of the irrotational part vanishes. In the 
irrotational (zero spin) case we have Di Uj — D, Ui = 



and thus ^^-'wi is derivable from a potential. For general 
rotation states we write 



i^)u'- ^D'(l) + w\ 



(33) 



so that DV and w^ denote the irrotational and rotational 
pieces of the velocity. In order to assure that w* is purely 
rotational one usually requires that 



Au;* = 0. 



(34) 



In subsection IIIFI we will show how one can choose w* 
such that Eq. (p4)) is satisfied. However, Eq. (p4)) is not 
explicitly used in any of the derivations in this subsection. 
Note that once ^'^'^Ui is known vP = —u^n^/a can be 
obtained from M^-Sp = —h?. If we choose w^n^ — the 
split of '^^fei in Eq. (|33p can be extended to 



r,^' — XJt^ 



V^(/. + w'', 



(35) 



where the time dependence of is now chosen such that 
it satisfies V°(/) = u" 

In order to simplify Eqs. (|30p and (pij) we now assume 
that 



£^{pou°) « £^9^,1, 







(36) 



but we will not assume that £(^Ui, vanishes as well. In- 
stead we assume that 



7r£e(v.^)«o, 



(37) 



so that the time derivative of the irrotational piece of the 
fluid velocity vanishes in corotating coordinates. Fur- 
thermore we also assume that 



^i£^w^ w 0, 
where we have defined 



e' = 



(38) 



(39) 



The assumption in Eq. (|38p describes the fact that the 
rotational piece of the fluid velocity (which gives rise to 
the spin) is constant along ^'^ which is parallel to the 
worldline of the star center. Defining 



A^M = ^A" _ ^M ^ (0^ A/fc* 



(40) 



and using Eqs. p7l) and p8l) the Lie derivative term in 
Eq. ((3T|) can be written as 



~ li£Ai'Wy^'^''^'>£/s,kWi. 



(41) 



Here ^'^'^£ is the Lie derivative in 3 dimensions. Thus 
Eqs. (150)) and ([5T|) simplify and can be rewritten as 



A (poau°l^ 



OtaA 







(42) 



and 



where we use the abbreviations 



In order to further simplify Eq. P5)) note that 






(43) 



(44) 



which follows from Eqs. (|29l), (|33l), ([39]) and (IMl). Hence 



'V+AkWi 



,-,0 



(3)£^u" + ^M3)£ ^^^^0 (45) 



where we have assumed that both vP and jik are ap- 
proximately constant along the 3-vector ^, which lies 
along the direction of the fluid's rotational velocity piece 
'u;^ Note, that ^'^^i^y+AfeWi is of order 0{w)^, while as- 
sumptions ([57]) and ([55]) are 0(1) and 0(w) in w*. Thus 
alternatively we can view Eq. (j45p as an assumption that 
will hold if w* is small compared to D^(j>. All three as- 
sumptions (|37)) . (|38t and (|45l) are discussed in appendix 
l^for a simple case. 

With the last assumption in Eq. (gS]) the Euler Eq. P5)) 
yields 



F^^'-Dfef/) = -C, 



(46) 



where C is a constant of integration, that is in general 
different for each star. 

In the corotating case where V^ = 0, Eq. (j42|) is iden- 
tically satisfied and Eq. pS)) reduces to 



h = -Cu", 

The u" here can be computed from u^u^ 
duces to 



u° - l/V«'-(/3.+?0(/3'+e) 

for V^ = 0. 

If the stars are not corotating V^ is given by 

In this case the continuity equation (|42l) becomes 



(47) 
-1 and re- 

(48) 



D, 



^^{D'cj) + w')-p^au\(i'+C) 



0. 



Note that u^u^ 



-1 yields 



u" = 



^K^ + {D,(t) + Wi){Di<j) + w') 



(49) 



(50) 



(51) 



so that Eq. ([50]) is a non-linear elliptic equation for 0. 
Using vP from Eq. ([5T|) the integrated Euler equation pS)) 
can then be solved for h with the result 



v/i2-(A</> + wO(i?^0 + u;*), 



(52) 



L' 



2^2 



i]2 



(53) 



and 



b = ae + l31D,<j> - Cf + 2a^{D,c^ + w^)w\ (54) 

Note that the rotational piece of the fluid velocity w* 
can be freely chosen, and that the fluid equations (|50|) 
and (15211 reduce to the well known result for irrotational 
stars [4l|, m if w^ = 0. 



F. Further simplifications and boundary conditions 

Next we also choose a maximal slice with K = Q, 
and assume that the conformal 3-metric is flat and given 
byiP 



lij 



(55) 



This latter assumption merely simplifies our equations 
and could in principle be improved by e.g. choosing 
a post-Newtonian expression for 7jj or by matching a 
post-Newtonian metric with a single neutron star solu- 
tion similar to |43l - l48| . Using Eq. (155]) the Hamiltonian 
and momentum constraints in Eq. © and Eq. ([TT)) sim- 
plify and we obtain 



D {aip) — aip 



32a- 



■(LB)*^(LB),j+27r^*(p + 25) 



(56) 



where (LBf^ = D'B^+D^B' - ^S^WkB^, Di = d^ and 



The elliptic equations (f56| have to be solved subject to 
the boundary conditions 



lim V' = 1, Inn B* — 0, lim aij) — 1 



(58) 



at spatial infinity. 

The equations ((56|) need to be solved together with the 
fluid equations (l5(I)l and ([5^ . These fluid equations sim- 
plify in corotating coordinates where ^' = 0. Further- 
more they can be expressed in terms of the derivative 
operator Di by noting that 



Di(t> = A0, D^(j) = V~ 


4l)V. 


(59) 


addition w* can be replaced by 






W' = TJJ-^W\ 




(60) 



The latter scaling is useful since 



(61) 



so that if we choose Diw'' = we automatically obtain 
Diw"^ = 0. One obvious choice for the conformal rota- 
tional velocity could be 



w 



ijk 



uj^x'' 



'-C* 



), 



(62) 



where x'q^ is the location of the star center, which could 
be defined as the point with the highest rest mass density 
Po or as the center of mass of the star. However, it is also 
possible to choose 



w 



/(k" 



^c* 



)e'^''LoHx'' 



^c*i^ 



(63) 



where /(|a;" — x^^\) is any function that only depends on 
the conformal distance from the star's center. Thus the 
method described here is capable if of giving an arbitrary 
rotational velocity to each star. 

Also note, that we need a boundary condition at the 
star surface to solve Eq. ([50]) . This boundary condi- 
tion can be obtained from Eq. ([5D|) itself by evaluating 
Eq. ([50l) on the boundary where po — >■ but Dipo ^ 0. 
Taking this limit we obtain 



{D'<j>)D,po + wW.po = hu°{P' + e)D,po 



(64) 



at the star surface. In applications it may be a good idea 
to choose w* such that w^Dipo vanishes, otherwise the 
rotational velocity has a component perpendicular to the 
star's surface. Also notice that Eq. ([50)1 together with its 
boundary condition in Eq. (j64p do not uniquely specify 
the solution. If solves both Eqs. ((5(1)) and ([M)) 0-1- const 
will be a solution as well. In numerical codes this kind of 
ambiguity is usually removed by adding e.g. the volume 
integral of 4> over the star to the boundary condition. 



III. THE NEWTONIAN LIMIT 

We now investigate Newtonian limit of the approxi- 
mate matter equations derived above. If ip is the New- 
tonian potential satisfying 9j9V = iwpo and u* — u^/u^ 
the Newtonian fluid velocity (in inertial coordinates) we 
can express the Newtonian limit as 



300 
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-l-2ip 
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-> 


l + ip 


9ot = A 


-> 





J^3 = l^l 


-^ 


S^, 
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-^ 


[n X xY 


i^^'^U, 


-^ 


Vi 


Vi 
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di<j) + Wi 
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-^ 


d'(h + w' 






= 50mW 
1 + hN 



-1- 



•^ 



P 

PO 



(65) 



where v — \/v^Vi and V^ is the fluid velocity in corotating 
coordinates. 

Using Eqs. (|65|) Eq. ^ reduces to 

d,{poV')^0 (66) 

which is the Newtonian continuity equation in corotating 
coordinates where dfPo — 0. 

In order to examine the limit of Eq. (|46p we first note 
that 



-f(3)^,l/fc 



-huf^e 



(67) 



D^{V''wk) = V'iDi^ - Dfuk) -^^^£vm. (68) 



and 



Using Eqs. ([S7)) and (p5| together with the limits in 
Eqs. (1551) the gradient of Eq. PS|) yields 

(69) 
In order to show that this is the Euler equation of New- 
tonian physics we first note that the time derivative dt' 
in corotating coordinates is related to the time derivative 
dt in inertial coordinates by df — dt -f'-'^'i^j. Then 

dt'Vi = df {di(j) + Wi- ^i) 

= df (^0) + dtW, +(3) £^w, - dt.^^ 

= df{d,(j)) ~ dt'^^ + (dtw, +(^)£^-w,) 

+^^^£v+AkW,~'^^'^£vw,. (70) 

In the last equality all terms but the last vanish, if we 
make the same assumptions as in Eqs. ([ST]) . ([55]) and 
(|^5)) . Hence Eq. ([S^ can be rewritten as 



dt'V^ + V''dkV^ + 2[nxVY + [nx {fix x)Y 



d,P 
Po 



^ 



(71) 

which is simply the well known Euler equation of Newto- 
nian physics expressed in corotating coordinates. Thus 
we see that our new matter equations reduce to the cor- 
rect result in the Newtonian limit. 



IV. DISCUSSION 

Realistic neutron stars in binaries will be spinning. 
From observations of millisecond pulsars we know that 
these spins can be substantial enough to influence the 
late inspiral and merger dynamics of the binary. 

There have been prior attempts to construct initial 
data for spinning neutron stars. In f21l| (hereafter MS) 



the Euler equation is not solved directly. Rather it is re- 
placed by an equation equivalent to ^ -t-^^^UfcT^*^ = —C. 
However, as already pointed out by MS, this equation 
agrees with the integrated Euler Eq. ([^S]) only for the 
corotating and the irrotational case. Thus in general the 
Euler equation is violated in the MS approach. Further- 
more, MS split u^/u^ and not ^'^fe* into an irrotational 
and a rotational part (see Eq. (j331))- This has two con- 
sequences. First, their equations do not have the correct 
limit in the irrotational case. And second, since u^ /vP is 
not a purely spatial vector it is inconsistent to set u^ /vP 
equal to something like D^cj) which is a purely spatial vec- 
tor. This explains why the continuity equation of MS has 
no shift terms unlike in Eq. ([50|) and in [4l|, H^] ■ When 
MS c omp are their results for a particular corotating case 
with [ 201 they find that their approach introduces errors 
of about 2% in the angular momentum. 

Another approach to include spin that is aligned with 
the orbital angular momentum was proposed in [26| . 
hereafter BS. This approach does not seek to analyti- 
cally integrate the Euler equation as we have done here. 
Instead the divergence of Eq. (|3ip is set to zero, which 
leads to another elliptic equation. However, as pointed 
out first by Gourgoulhon ^4^ , in general the Euler equa- 
tion itself is not satisfied if we only enforce its divergence 
to be zero. Hence the BS approach can lead to initial 
data that do not obey the Euler equation. Furthermore, 
the boundary condition given by BS for their new elliptic 
equation seems to imply that the star surface is always 
at the same location. If we consider the usual numerical 
treatment where we start from an initial guess for the 
stars which is iteratively refined, it is unclear how the 
star surface can change during the iterations. 

The purpose of this paper is thus to introduce a new 
method for the computation of binary neutron star initial 
data with arbitrary rotation states. Our method is de- 
rived from the standard matter equations of perfect fluids 
together with certain quasi-equilibrium assumptions. We 
assume that there is an approximate helical Killing vector 
^^ and that Lie derivatives of the metric variables with 
respect to ^^ vanish. We also assume that scalar matter 
variables such as h or po have Lie derivatives that vanish 
with respect to ^^ . However, as discussed in appendix 
1X1 the Lie derivative of the fluid velocity u^ is expected 
to be non-zero for arbitrary spins. We split the fluid 
velocity u^ into an irrotational piece (derived from a po- 
tential 0) and a rotational piece w^ , and assume that only 
the irrotational piece has a vanishing Lie derivative (see 
Eq. (1571) ) with respect to ^^ . This can be interpreted the 
natural generalization of the irrotational case where one 
commonly assumes £^hu'^ = 0. Furthermore we know 
that the spin of each star remains approximately constant 
since the viscosity of the stars is insufficient for tidal cou- 
pling [231 • To incorporate this fact, we use Eq. ([551) which 
is based on the assumption that w* is constant along the 
star's motion described by the irrotational velocity piece 
Vc/). Since V^t/) is equivalent to the velocity of the star 
center, this latter assumption captures the fact that the 



spin or rotational velocity w'^ of each star remains ap- 
proximately constant. With these two assumptions the 
Euler equation simplifies to Eq. ([43| . In order to analyti- 
cally integrate Eq. ([^5|) we use the additional assumption 
([45[) that vP and 7y are constant along the field lines of 
the rotational velocity piece. We then arrive at the two 
matter equations Eq. ([SHI) and ([5^ . These equations 
reduce to well known equations 41, 42] for the irrota- 
tional case of Wi = 0. They also reduce to the corotating 
limit (where V* = 0) as is evident from Eqs. ([^^ and 
(|46p which are written in terms of V^. Furthermore, our 
equations reduce to the correct Newtonian limit. 

The elliptic equation in Eq. (|50|) can to solved (for </>) 
together with the Eqs. ([551) fo'" the metric variables once 
the enthalpy h is known. However, the enthalpy given 
by Eq. (|52p depends on the metric variables, and w^. 
Apart from their dependence on w^ this set of equations 
has a similar structure as for the case of irrotational neu- 
tron stars (where w' vanishes). The standard way (see 
e.g. |lfl]) to solve such a mixture of elliptic and alge- 
braic equations is by iteration, where at each step we 
first solve the elliptic equations for a given h and then 
use the algebraic Eq. (|5^ to update h. At each step we 
also need to specify w^. One way to do this would be by 
choosing a constant iD* as in Eq. ([5^. Note however, that 
other choices for w* are possible. We plan to investigate 
these possibilities in future numerical studies of our new 
method. For such studies it might useful to use a numeri- 
cal code like LORENE [ll,|5il3 or SGRID fll,[5l|53, 
where the star surface is always at a domain boundary 
so that the boundary condition in Eq. (IMl) can be easily 
implemented. 
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Appendix A: The matter quasi-equilibrium 
assumptions in a simplified case 

In the Newtonian limit ^■^•'{Ji is the equal to the fluid 
velocity in the inertial frame. When the two stars are well 
separated it is clear that each star is well approximated 
by an orbiting and spinning sphere. In this case the fluid 
velocity inside a star is given by 



'-^\ Ki [Q X xc*]i + [uj X (x ~ xc*)]i, 



(Al) 



where xc*i is the (time dependent) location of the star 
center, and Vli and Ui are the angular velocities of the 
orbital and spinning motion of the star. Within this ap- 
proximation we then have 



(j)^ [fl X Xc*]kx'' 



Wi 



[uj X {x ~ xc*)] 



(A2) 
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It is then easy to verify that the assumptions in Eqs. ((37)) 
and ([55]) are identically satisfied. Furthermore for ap- 
proximate spherical symmetry we see that the assump- 
tions in Eq. (|45p hold as well. In addition, we find that 

Afci ^[nx x], - D,(j) = [fi X (a; - xc*%. (A3) 

From this it follows that 

^X^iUy ='^^^ £i\h'Wi = (r^iWj - iOiVLj){x^ - x^fjj (A4) 



which illustrates that £^Ui, does not vanish even in this 
simplified case. The only case when £^Ui, can vanish is if 
the spin is aligned with the orbital angular momentum, 
i.e. if ujj = aQj for some constant a. 
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